#include "header.h"

static int N, n2;
static FILE *fh_out;
static __float128 *k, g;
static fftwq_complex *a[4], *RS;
static aux_ptr ax;

void init_prep(params_ptr in) {
  N = in->n; n2 = N/2 + 1; 
  //u = fftwq_malloc(sizeof(__float128)*N);
  k = fftwq_malloc(sizeof(__float128)*N);
  for (int j = 0; j < in->n; j ++) {
    //u[j] = 2.Q*pi*j/N;
    k[j] = j;
    if (j > N/2) k[j] = j - N;
  }
  for (int j = 0; j < 4; j++) 
  g = in->g;
}

void init_full_write(params_ptr in, aux_ptr in2) {
  ax = in2;
  ax->XX = fftwq_malloc(sizeof(__float128)*(N+1));
  ax->YY = fftwq_malloc(sizeof(__float128)*(N+1));
}

void write_out(fftwq_complex* in, char* str) {
  fh_out = fopen(str, "w");
  fprintf(fh_out, "#1.u 2.k 3.re R 4. im R\n\n");
  for (int j = 0; j < N; j++) {
    fprintf(fh_out, "%.15Qe\t%.15Qe\t%.15Qe\t%.15Qe\n", 2.Q*pi*j/N, k[j], 
                        creal(in[j]), cimag(in[j]) );
  }
  fclose(fh_out);
}

void writer_out(__float128* in, char* str) {
  fh_out = fopen(str, "w");
  fprintf(fh_out, "#1.u 2.k 3.re R 4. im R\n\n");
  for (int j = 0; j < n2; j++) {
    fprintf(fh_out, "%.28Qe\t%.28Qe\t%.28Qe\t%.28Qe\n", 2.Q*pi*j/N, k[j], 0., in[j] );
  }
  for (int j = 0; j < n2-2; j++) fprintf(fh_out, "%.28Qe\t%.28Qe\t%.28Qe\t%.28Qe\n", 2.Q*pi*(j+n2)/N, k[j+n2], 0.Q ,in[n2-2-j]);
  fclose(fh_out);
}

void writer_full(work_ptr wrk, char* str) {
  fh_out = fopen(str, "w");
  fprintf(fh_out, "#1.u 2.k 3. x-tilde 4. y 5. yk \n\n");
  fprintf(fh_out, "%.28Qe\t%Qe\t%.28Qe\t%.28Qe\t%.28Qe\n", 0.Q, 0.Q, 0.Q, wrk->YY[0], wrk->YYk[0]);
  for (int j = 1; j < n2-1; j++) {
    fprintf(fh_out, "%.28Qe\t%Qe\t%.28Qe\t%.28Qe\t%.28Qe\n", 2.Q*pi*j/N, k[j], wrk->XX[j-1], wrk->YY[j], wrk->YYk[j]);
  } 
  fprintf(fh_out, "%.28Qe\t%Qe\t%.28Qe\t%.28Qe\t%.28Qe\n", pi, k[n2-1], 0.Q, wrk->YY[n2-1], wrk->YYk[n2-1]);
  // must write other half
  for (int j = 1; j < n2-1; j++) fprintf(fh_out, "%.28Qe\t%Qe\t%.28Qe\t%.28Qe\t%.28Qe\n", 2.Q*pi*(j+n2-1)/N, k[j+n2-1], -1.Q*(wrk->XX[n2-2-j]),
      wrk->YY[n2-1-j], wrk->YYk[n2-1-j]);
  fclose(fh_out);
}

void write_ImS(work_ptr wrk) {
  fh_out = fopen("rawImS_all.dat","w");
  fprintf(fh_out, "\n");
  for (int j = 0; j < N+1; j++) fprintf(fh_out, "%.31Qe\n", 2.Q*pi*j/N - pi);
  fclose(fh_out);
}

void write_Rew(work_ptr wrk) {
  //fh_out = fopen("rawRew_all.dat","w");  
  __float128 avg = wrk->YY[0];
  //fprintf(fh_out, "\n");
  for (int j = 0; j < N+1; j++) {
    if (j < n2) {
       if (j == 0) {
         //fprintf(fh_out, "%.31Qe\n", 0.Q);
	 ax->YY[j] = 0.Q;
       }
       if (j != 0) {
         //fprintf(fh_out, "%.31Qe\n", avg - 1.Q*wrk->YY[j]);
	 ax->YY[j] = avg - 1.Q*wrk->YY[j];
       }
    }
    if (j >= n2) {
	if (j != N) {
	  //fprintf(fh_out, "%.31Qe\n", avg - 1.Q*wrk->YY[N - j]);
	  ax->YY[j] = avg - 1.Q*wrk->YY[N - j];
	}
        if (j == N) {
	  //fprintf(fh_out, "%.31Qe\n", 0.Q );
	  ax->YY[j] = 0.Q;
	}
    }
  }
  //fclose(fh_out);
}

void write_Imw(work_ptr wrk) {
  //fh_out = fopen("rawImw_all.dat","w");
  //fprintf(fh_out, "\n");
  for (int j = 0; j < N+1; j++) {
    if (j < n2-1) {
      if (j == 0){
        //fprintf(fh_out, "%.31Qe\n", 0.Q);
	ax->XX[j] = 0;
      }
      if (j != 0) {
	//fprintf(fh_out, "%.31Qe\n", 1.Q*wrk->XX[j-1]);
	ax->XX[j] = 1.Q*wrk->XX[j-1];
      }
    }
    if (j == n2-1) {
	//fprintf(fh_out, "%.31Qe\n", 0.Q);
	ax->XX[j] = 0.Q;	
    }
    if (j > n2-1) {
	if (j != N) {
	  //fprintf(fh_out, "%.31Qe\n", -1.Q*wrk->XX[N - j - 1]);
	  ax->XX[j] = -1.Q*wrk->XX[N - j - 1];	
	}
        if (j == N) {
	  //fprintf(fh_out, "%.31Qe\n", 0.Q );
	  ax->XX[j] = 0.Q;	
	}
    }
  }
  //fclose(fh_out);
}



